Pylops - Marchenko redatuming¶
Author: M.Ravasi¶
In this notebook we consider to extensions of the Marchenko inversion in the case of band-limitied reflectivity (i.e., convolved with wavelet).
We modify the Marchenko equations to account for it and we also use a broad-band (i.e., spiky) initial focusing function.
$$ \begin{bmatrix} \Theta \mathbf{P} \mathbf{f_d^+} \\ \mathbf{0} \end{bmatrix} = \begin{bmatrix} \mathbf{W} & \Theta \mathbf{P} \\ \Theta \mathbf{P^*} & \mathbf{W^*} \end{bmatrix} \begin{bmatrix} \mathbf{f^-} \\ \mathbf{f_m^+} \end{bmatrix} $$We also consider now the alternating problem of estimating both the focusing functions and the wavelet, using the following equation to update the wavelet.
$$ \begin{bmatrix} \Theta \mathbf{P} \mathbf{f^+} \\ \Theta \mathbf{P^*} \mathbf{f^-} \end{bmatrix} = \begin{bmatrix} \mathbf{F^-} \\ \mathbf{F_m^+} \mathbf{S} \end{bmatrix} \mathbf{w} \rightarrow \mathbf{d_w} = \mathbf{F_w} \mathbf{w} $$where $\mathbf{P}$ is the bandlimited reflection response (i.e. $P(t, r, s) = w(t) * P(t, r, s)$), $\mathbf{F^-}$ and $\mathbf{F_m^+}$ are convolution operators and $\mathbf{S}$ is a flipping operator that flips the time axis.
In [1]:
%load_ext memory_profiler
%load_ext autoreload
%autoreload 2
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.sparse import csr_matrix, vstack
from scipy.linalg import lstsq, solve
from scipy.sparse.linalg import cg, lsqr
from scipy.signal import convolve, fftconvolve, filtfilt, hilbert
from pylops import LinearOperator
from pylops.utils import dottest
from pylops.utils.wavelets import *
from pylops.utils.seismicevents import *
from pylops.utils.tapers import *
from pylops.basicoperators import *
from pylops.signalprocessing import *
from pylops.waveeqprocessing.mdd import *
from pylops.waveeqprocessing.marchenko import *
from pylops.optimization.leastsquares import *
from pylops.optimization.sparsity import *
Inputs¶
Input parameters
In [2]:
inputfile = '../data/marchenko/input.npz' # choose file in testdata folder of repo
vel = 2400.0 # velocity
toff = 0.045 # direct arrival time shift
nsmooth = 10 # time window smoothing
nfmax = 500 # max frequency for MDC (#samples)
nstaper = 11 # source/receiver taper lenght
n_iter = 20 # iterations
jr = 1 # subsampling in r
js = 1 # subsampling in s
clip = 5e5 # figure clipping
Load input
In [3]:
inputdata = np.load(inputfile)
Read and visualize geometry
In [4]:
# Receivers
r = inputdata['r'][:,::jr]
nr = r.shape[1]
dr = r[0,1]-r[0,0]
# Sources
s = inputdata['s'][:,::js]
ns = s.shape[1]
ds = s[0,1]-s[0,0]
# Virtual points
vs = inputdata['vs']
# Density model
rho = inputdata['rho']
z, x = inputdata['z'], inputdata['x']
plt.figure(figsize=(18,9))
plt.imshow(rho, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker='*', s=150, c='r', edgecolors='k')
plt.scatter(r[0, ::10], r[1, ::10], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(vs[0], vs[1], marker='.', s=250, c='m', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]'),plt.title('Model and Geometry')
plt.xlim(x[0], x[-1]);
Read data
In [5]:
# time axis
t = inputdata['t']
t2 = np.concatenate([-t[::-1], t[1:]])
ot, dt, nt, nt2 = t[0], t[1]-t[0], len(t), len(t2)
# data
R = inputdata['R'][::js, ::jr]
R = np.swapaxes(R, 0, 1) # R[s, r, f] assume particle velocity receiver and integrate over them (via reciprocity)
# tapering
taper = taper3d(nt, [ns, nr], [nstaper, nstaper], tapertype='hanning')
R = R*taper
Define wavelet
In [6]:
# wavelet
ntwav = 21
f0 = 20
phi = 0
amp = .7
ampf = 0.
# original (simple ricker)
wavsymm, twav, wav_c = ricker(t[:ntwav], f0)
# modified wavelet
wav = amp * wavsymm
WAV = np.fft.rfft(wav, 2**10)
WAVabs = np.abs(WAV)
WAVabsmod = WAVabs + filtfilt(np.ones(10)/10., 1, np.random.normal(0., ampf, len(WAV)), method='pad')
WAV = WAVabsmod * np.exp(1j*np.angle(WAV))
wav = np.fft.irfft(WAV, 2**10)[:len(wav)]
if phi != 0:
wav = amp * wavsymm(np.cos(np.deg2rad(phi))*wavsymm + np.sin(np.deg2rad(phi))*np.imag(hilbert(wavsymm)))
plt.figure(figsize=(8, 3))
plt.plot(WAVabs, 'k')
plt.plot(WAVabsmod, '--r')
plt.title('Wavelet');
plt.figure(figsize=(8, 3))
plt.plot(wavsymm, 'k', label='Original')
plt.plot(wav, 'r', label='Modified')
plt.legend()
plt.title('Wavelet');
Apply wavelet to R
In [7]:
# Convolve R with wavelet
#Rwav = np.apply_along_axis(convolve, 2, R, wav, mode='full')
#Rwav = Rwav[:,:,wav_c:][:, :, :nt]
Rwavsymm = fftconvolve(R, wavsymm[np.newaxis,np.newaxis,:], mode='same', axes=-1)
Rwav = fftconvolve(R, wav[np.newaxis,np.newaxis,:], mode='same', axes=-1)
fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(R[ns//2].T, cmap='gray', vmin=-1e-2, vmax=1e-2, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('R shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
ax.set_ylim(1.5, 0)
fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(Rwavsymm[ns//2].T, cmap='gray', vmin=-1e-1, vmax=1e-1, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('Rwavsymm shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
ax.set_ylim(1.5, 0)
fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(Rwav[ns//2].T, cmap='gray', vmin=-1e-1, vmax=1e-1, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('Rwav shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight');
ax.set_ylim(1.5, 0);
Read subsurface fields
In [8]:
Gsub = inputdata['Gsub'][:, ::jr]
G0sub = inputdata['G0sub'][:, ::jr]
# convolve G with wavelet
Gsub = np.apply_along_axis(convolve, 0, Gsub, wav, mode='full')
Gsub = Gsub[wav_c:][:nt]
In [9]:
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub/Gsub.max(), cmap='gray', vmin=-1, vmax=1,
extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax2.imshow(G0sub/G0sub.max(), cmap='gray', vmin=-1, vmax=1,
extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax2.set_title(r'$G0$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(G0sub[:, nr//2]/G0sub.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
In [10]:
# Direct arrival window - traveltime
toffw = toff - 0.015
directVS = np.sqrt((vs[0]-r[0])**2+(vs[1]-r[1])**2)/vel
directVS_off = directVS - toffw
# Window
idirectVS_off = np.round(directVS_off/dt).astype(np.int)
w = np.zeros((nr, nt))
for ir in range(nr):
w[ir, :idirectVS_off[ir]]=1
w = np.hstack((np.fliplr(w), w[:, 1:]))
if nsmooth>0:
smooth=np.ones(nsmooth)/nsmooth
w = filtfilt(smooth, 1, w)
fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
im = ax.imshow(w.T, cmap='gray', extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax.plot(r[0], directVS,'b'),ax.plot(r[0], -directVS,'b')
ax.set_title('Window'), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
fig.colorbar(im, ax=ax);
Inversion of Marchenko equations¶
Use same wavelet for Wop and R¶
Inversion with modified wavelet
In [11]:
# Add negative time to operators
Rtwosided = np.concatenate((np.zeros((ns, nr, nt-1)), Rwav), axis=-1)
Rtwosided_fft = np.fft.rfft(Rtwosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
Rtwosided_fft = Rtwosided_fft[...,:nfmax]
Rtwosided_fft = Rtwosided_fft.transpose(2, 0, 1) # R[f, s, r]
# Operators
Rop = MDC(Rtwosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=dr,
twosided=True, transpose=False, dtype='complex64')
R1op = MDC(Rtwosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=dr,
twosided=True, transpose=False, conj=True, dtype='complex64')
Rollop = Roll((2*nt-1) * ns,
dims=(2*nt-1, ns),
dir=0, shift=-1)
Wavop = Convolve1D(ns*nt2, wav, wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Wav1op = Convolve1D(ns*nt2, np.flipud(wav), wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Wop = Diagonal(w.T.flatten())
# Input focusing function
fd_plus = np.concatenate((np.fliplr(G0sub.T).T, np.zeros((nt-1, nr))))
In [12]:
Mop = VStack([HStack([Wavop, -1*Wop*Rop]),
HStack([-1*Wop*Rollop*R1op, Wav1op])])*BlockDiag([Wop, Wop])
Gop = VStack([HStack([Wavop, -1*Rop]),
HStack([-1*Rollop*R1op, Wav1op])])
dottest(Gop, 2*ns*(2*nt-1), 2*nr*(2*nt-1), verb=True)
dottest(Mop, 2*ns*(2*nt-1), 2*nr*(2*nt-1), verb=True);
In [13]:
p0_minus = Rop * fd_plus.flatten()
p0_minus = p0_minus.reshape(nt2, ns).T
d = Wop*Rop*fd_plus.flatten()
d = np.concatenate((d.reshape(nt2, ns), np.zeros((nt2, ns))))
f1_adj = Mop.H*d.flatten()
f1_inv = lsqr(Mop, d.ravel(), iter_lim=n_iter*2, show=True)[0]
res = d.ravel() - Mop * f1_inv
f1_adj = f1_adj.reshape(2*nt2, nr)
f1_inv = f1_inv.reshape(2*nt2, nr)
res = res.reshape(2*nt2, nr)
tau = np.linalg.norm(f1_inv.ravel(), 1)
print('|f1_inv|_1=', tau)
f1_adj_tot = f1_adj + np.concatenate((np.zeros((nt2, nr)), fd_plus))
f1_inv_tot = f1_inv + np.concatenate((np.zeros((nt2, nr)), fd_plus))
In [14]:
g_adj = Gop*f1_adj_tot.flatten()
g_inv = Gop*f1_inv_tot.flatten()
g_adj = g_adj.reshape(2*nt2, ns)
g_inv = g_inv.reshape(2*nt2, ns)
f1_adj_minus, f1_adj_plus = f1_adj_tot[:nt2].T, f1_adj_tot[nt2:].T
f1_inv_minus, f1_inv_plus = f1_inv_tot[:nt2].T, f1_inv_tot[nt2:].T
g_adj_minus, g_adj_plus = -g_adj[:nt2].T, np.fliplr(g_adj[nt2:].T)
g_inv_minus, g_inv_plus = -g_inv[:nt2].T, np.fliplr(g_inv[nt2:].T)
g_inv_tot = g_inv_minus + g_inv_plus
In [15]:
# Need to recreate combined data as new implementation stacks over time instead of space
d_plot = np.concatenate((d[:(2*nt-1)], d[(2*nt-1):]), axis=1).T
res_plot = np.concatenate((res[:(2*nt-1)], res[(2*nt-1):]), axis=1).T
f1_adj_tot_plot = np.concatenate((f1_adj_tot[:(2*nt-1)], f1_adj_tot[(2*nt-1):]), axis=1).T
f1_inv_tot_plot = np.concatenate((f1_inv_tot[:(2*nt-1)], f1_inv_tot[(2*nt-1):]), axis=1).T
g_adj_plot = np.concatenate((g_adj[:(2*nt-1)], g_adj[(2*nt-1):]), axis=1).T
g_inv_plot = np.concatenate((g_inv[:(2*nt-1)], g_inv[(2*nt-1):]), axis=1).T
fig, axs = plt.subplots(1, 6, sharey=True, figsize=(16, 7))
axs[0].imshow(d_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('d'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1)
axs[1].imshow(res_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title('res'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1)
axs[2].imshow(f1_adj_tot_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(f1_inv_tot_plot.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[3].set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[3].axis('tight')
axs[3].set_ylim(1, -1);
axs[4].imshow(g_adj_plot.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[4].set_title(r'$g_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[4].axis('tight')
axs[4].set_ylim(1, -1);
axs[5].imshow(g_inv_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[5].set_title(r'$g_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[5].axis('tight')
axs[5].set_ylim(1, -1);
In [16]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_adj_minus.T, cmap='gray', vmin=-10*clip, vmax=10*clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_adj_plus.T, cmap='gray', vmin=-10*clip, vmax=10*clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [17]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
Use different wavelet for Wop and R¶
Let's what happens if we use the wrong wavelet (use wavsymm for the wavelet operator but R with the modified wavelet)
In [18]:
Wavsymmop = Convolve1D(ns*nt2, wavsymm, wav_c, dims=(nt2, ns),
dir=0, dtype='float64')
Wavsymm1op = Convolve1D(ns*nt2, np.flipud(wavsymm), wav_c, dims=(nt2, ns),
dir=0, dtype='float64')
# Input focusing function
fd_plus = np.concatenate((np.fliplr(G0sub.T).T, np.zeros((nt-1, nr))))
Msymmop = VStack([HStack([Wavsymmop, -1*Wop*Rop]),
HStack([-1*Wop*Rollop*R1op, Wavsymm1op])])*BlockDiag([Wop, Wop])
Gsymmop = VStack([HStack([Wavsymmop, -1*Rop]),
HStack([-1*Rollop*R1op, Wavsymm1op])])
p0_symm_minus = Rop * fd_plus.flatten()
p0_symm_minus = p0_symm_minus.reshape(nt2, ns).T
d_symm = Wop * Rop * fd_plus.flatten()
d_symm = np.concatenate((d_symm.reshape(nt2, ns), np.zeros((nt2, ns))))
f1_symm_adj = Msymmop.H*d_symm.flatten()
f1_symm_inv = lsqr(Msymmop, d_symm.ravel(), iter_lim=n_iter*2, show=True)[0]
res_symm = d.ravel() - Msymmop * f1_symm_inv
f1_symm_adj = f1_symm_adj.reshape(2*nt2, nr)
f1_symm_inv = f1_symm_inv.reshape(2*nt2, nr)
res_symm = res_symm.reshape(2*nt2, nr)
f1_symm_adj_tot = f1_symm_adj + np.concatenate((np.zeros((nt2, nr)), fd_plus))
f1_symm_inv_tot = f1_symm_inv + np.concatenate((np.zeros((nt2, nr)), fd_plus))
g_symm_adj = Gsymmop*f1_symm_adj_tot.flatten()
g_symm_inv = Gsymmop*f1_symm_inv_tot.flatten()
g_symm_adj = g_symm_adj.reshape(2*nt2, ns)
g_symm_inv = g_symm_inv.reshape(2*nt2, ns)
f1_symm_adj_minus, f1_symm_adj_plus = f1_symm_adj_tot[:nt2].T, f1_symm_adj_tot[nt2:].T
f1_symm_inv_minus, f1_symm_inv_plus = f1_symm_inv_tot[:nt2].T, f1_symm_inv_tot[nt2:].T
g_symm_adj_minus, g_symm_adj_plus = -g_symm_adj[:nt2].T, np.fliplr(g_symm_adj[nt2:].T)
g_symm_inv_minus, g_symm_inv_plus = -g_symm_inv[:nt2].T, np.fliplr(g_symm_inv[nt2:].T)
g_symm_inv_tot = g_symm_inv_minus + g_symm_inv_plus
In [19]:
# Need to recreate combined data as new implementation stacks over time instead of space
res_plot = np.concatenate((res_symm[:(2*nt-1)], res_symm[(2*nt-1):]), axis=1).T
f1_adj_tot_plot = np.concatenate((f1_symm_adj_tot[:(2*nt-1)], f1_symm_adj_tot[(2*nt-1):]), axis=1).T
f1_inv_tot_plot = np.concatenate((f1_symm_inv_tot[:(2*nt-1)], f1_symm_inv_tot[(2*nt-1):]), axis=1).T
g_adj_plot = np.concatenate((g_symm_adj[:(2*nt-1)], g_symm_adj[(2*nt-1):]), axis=1).T
g_inv_plot = np.concatenate((g_symm_inv[:(2*nt-1)], g_symm_inv[(2*nt-1):]), axis=1).T
fig, axs = plt.subplots(1, 6, sharey=True, figsize=(16, 7))
axs[0].imshow(d_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('d'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1)
axs[1].imshow(res_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title('res'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1)
axs[2].imshow(f1_adj_tot_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(f1_inv_tot_plot.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[3].set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[3].axis('tight')
axs[3].set_ylim(1, -1);
axs[4].imshow(g_adj_plot.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[4].set_title(r'$g_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[4].axis('tight')
axs[4].set_ylim(1, -1);
axs[5].imshow(g_inv_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[5].set_title(r'$g_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[5].axis('tight')
axs[5].set_ylim(1, -1);
In [20]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_symm_adj_minus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_symm_adj_plus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_symm_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_symm_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [21]:
fig, ax = plt.subplots(1, 1, figsize=(12, 2))
ax.plot(Rwavsymm[50, 50], 'k', label='Original')
ax.plot(Rwav[50, 50], 'r', label='Modified')
ax.set_title('R')
ax.set_xlim(0, 300)
ax.legend()
fig, axs = plt.subplots(2, 1, figsize=(12, 5))
axs[0].plot(f1_symm_inv_minus[50], 'k')
axs[0].plot(f1_inv_minus[50], 'r')
axs[0].set_title('f^-')
axs[0].set_xlim(600, 900)
axs[1].plot(f1_symm_inv_plus[50], 'k', label='Original')
axs[1].plot(f1_inv_plus[50], 'r', label='Modified')
axs[1].set_title('f^+')
axs[1].set_xlim(600, 900)
axs[1].legend();
In [22]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_symm_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_symm_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_symm_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_symm_inv_tot[nr//2, nt-1:]/g_symm_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
We compare now the results
In [23]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(f1_adj_minus.T - f1_symm_inv_minus.T, cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$f^-$ error adj'), axs[0].set_xlabel(r'$x_R$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_adj_plus.T - f1_symm_adj_plus.T, cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^+$ error adj'), axs[1].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
In [24]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(f1_inv_minus.T - f1_symm_inv_minus.T, cmap='gray',
vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$f^-$ error inv'), axs[0].set_xlabel(r'$x_R$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_plus.T - f1_symm_inv_plus.T, cmap='gray',
vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^+$ error inv'), axs[1].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
In [25]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(g_inv_minus.T - g_symm_inv_minus.T, cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$g^-$ error inv'), axs[0].set_xlabel(r'$x_R$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_plus.T + g_symm_inv_plus.T, cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^+$ error inv'), axs[1].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
We can see how in both cases we have been able to come up with a set of focusing functions that can explain the data (near to zero residual).
However only the first solution is correct and leads to no artefacts in $G^-$ while the second does. Althought not always valid, we can perhaps try to use this argument (minimum L2 norm of $G^-$) as a regularization term when inverting for the wavelet?
Wavelet estimation equations¶
Let's try to invert for the wavelet given the already estimated focusing functions
In [26]:
f1_inv_tot = f1_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_inv_minus, f1_inv_plus = f1_inv_tot[:nt2].T, f1_inv_tot[nt2:].T
d1 = Wop*Rop*f1_inv_plus.T.flatten()
d2 = Wop*R1op*f1_inv_minus.T.flatten()
dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))
In [27]:
W1op = Diagonal(w.flatten())
F1plus = VStack([Convolve1D(nt2, f1_inv[nt2:, ir], offset=nt-2) for ir in range(nr)])
F1minus = VStack([Convolve1D(nt2, f1_inv[:nt2, ir], offset=nt-1) for ir in range(nr)])
S = Flip(nt2)
Fw = VStack([W1op*F1minus, W1op*F1plus*S])
In [28]:
extended_wav = np.pad(wav, (nt-ntwav, nt-ntwav), mode='constant')
dwpred = Fw * extended_wav
dwpred = dwpred.reshape(2 * nr, nt2)
In [29]:
fig, axs = plt.subplots(1, 3, figsize=(15, 9))
axs[0].imshow(dw.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[0].axis('tight');
axs[1].imshow(dwpred.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[1].axis('tight')
axs[2].imshow(dw.T - dwpred.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[2].axis('tight');
plt.figure(figsize=(12, 3))
plt.plot(dw[50])
plt.plot(dwpred[50]);
plt.xlim(600, 800);
plt.figure(figsize=(12, 3))
plt.plot(dw[150])
plt.plot(dwpred[150]);
plt.xlim(600, 800);
In [30]:
wavest = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0]
In [31]:
plt.figure()
plt.plot(twav, wavsymm, '--k', lw=1, label='Modified (used)')
plt.plot(twav, extended_wav[nt-ntwav:nt+ntwav-1], 'k', lw=2, label='Original')
plt.plot(twav, wavest[nt-ntwav:nt+ntwav-1], '--r', label='Estimated')
plt.grid()
You can see that the correct wavelet has been estimated.
We now use the focusing functions retrieved using the the wrong wavelet (which could have been the initial guess in the alternating algorithm)
In [32]:
f1_symm_inv_tot = f1_symm_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_symm_inv_minus, f1_symm_inv_plus = f1_symm_inv_tot[:nt2].T, f1_symm_inv_tot[nt2:].T
d1_symm = Wop*Rop*f1_symm_inv_plus.T.flatten()
d2_symm = Wop*R1op*f1_symm_inv_minus.T.flatten()
dw_symm = np.concatenate((d1_symm.reshape(nt2, ns).T, d2_symm.reshape(nt2, ns).T))
W1op = Diagonal(w.flatten())
F1plus = VStack([Convolve1D(nt2, f1_symm_inv[nt2:, ir], offset=nt-2) for ir in range(nr)])
F1minus = VStack([Convolve1D(nt2, f1_symm_inv[:nt2, ir], offset=nt-1) for ir in range(nr)])
S = Flip(nt2)
Fw = VStack([W1op*F1minus, W1op*F1plus*S])
extended_wav = np.pad(wavsymm, (nt-ntwav, nt-ntwav), mode='constant')
dwpred_symm = Fw * extended_wav
dwpred_symm = dwpred_symm.reshape(2 * nr, nt2)
In [33]:
fig, axs = plt.subplots(1, 3, figsize=(15, 9))
axs[0].imshow(dw_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[0].axis('tight');
axs[1].imshow(dwpred_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[1].axis('tight')
axs[2].imshow(dw_symm.T - dwpred_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[2].axis('tight');
plt.figure(figsize=(12, 3))
plt.plot(dw_symm[50])
plt.plot(dwpred_symm[50]);
plt.xlim(600, 800);
plt.figure(figsize=(12, 3))
plt.plot(dw_symm[150])
plt.plot(dwpred_symm[150]);
plt.xlim(600, 800);
In [34]:
wavest_symm = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0]
plt.figure()
plt.plot(twav, wavsymm, '--k', lw=1, label='Modified (used)')
plt.plot(twav, extended_wav[nt-ntwav:nt+ntwav-1], 'k', lw=2, label='Original')
plt.plot(twav, wavest_symm[nt-ntwav:nt+ntwav-1], '--r', label='Estimaged')
plt.grid()
We can see that the inversion returns a completely new wavelet when using the wrong focusing functions. This means that if we let the focusing functions converge (zero residual) and estimate the wavelet at that point we will get out something that is compensating for the fact that we have used incosistent wavelets in the data and Wop.
We need to try using a sparse solver for the focusing functions to have slow convergence and use an alternating optimization as in EPSI.
Here we just try the sparse solver and let it converge
In [52]:
# Initial focusing function
iG0 = np.argmax(G0sub, axis=0)
G0sub_sparse = np.zeros_like(G0sub)
G0sub_sparse[iG0, np.arange(nr)] = G0sub.max()
G0sub_sparse = fftconvolve(G0sub_sparse, np.hanning(7)[:, np.newaxis], mode='same', axes=0)
#G0sub_sparse = G0sub.copy()
fd_plus_sparse = np.concatenate((np.fliplr(G0sub_sparse.T).T, np.zeros((nt-1, nr))))
d_sparse = Wop*Rop*fd_plus_sparse.flatten()
d_sparse = np.concatenate((d_sparse.reshape(nt2, ns), np.zeros((nt2, ns))))
In [53]:
f1_sparse_inv, _, info = SPGL1(Mop, d_sparse.flatten(), tau=tau, iter_lim=100, verbosity=2)
f1_sparse_inv = f1_sparse_inv.reshape(2*nt2, nr)
print('|f1_inv|_1=', np.linalg.norm(f1_sparse_inv.ravel(), 1))
f1_sparse_inv_tot = f1_sparse_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus_sparse))
In [54]:
plt.figure()
plt.plot(info['xnorm1'])
plt.axhline(tau, color='r');
In [55]:
g_sparse_inv = Gop*f1_sparse_inv_tot.flatten()
g_sparse_inv = g_sparse_inv.reshape(2*nt2, ns)
f1_sparse_inv_minus, f1_sparse_inv_plus = f1_sparse_inv_tot[:nt2].T, f1_sparse_inv_tot[nt2:].T
g_sparse_inv_minus, g_sparse_inv_plus = -g_sparse_inv[:nt2].T, np.fliplr(g_sparse_inv[nt2:].T)
g_sparse_inv_tot = g_sparse_inv_minus + g_sparse_inv_plus
In [56]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_sparse_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_sparse_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [57]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_sparse_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_sparse_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_sparse_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_sparse_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
Alternating optimization¶
Finally we perform alternating optimization to estimate both focusing functions and wavelet
In [64]:
optimalsigma = 0.01 * np.linalg.norm(d.ravel())
nouter = 10
quickplot=True
# initial wavelet
ntwav = len(twav)//2 + 1
nfft = 2**10
In [65]:
# Estimate wavelet spectrum
wavest0_fft = np.mean(np.abs(np.fft.fft(Rwav, nfft, axis=-1)), axis=(0, 1))
fwest = np.fft.fftfreq(nfft, d=dt/1000)
# Create wavelet in time
wavest0 = np.real(np.fft.ifft(wavest0_fft)[:ntwav])
wavest0 = np.concatenate((np.flipud(wavest0[1:]), wavest0), axis=0)
wavest0 = wavest0 / wavest0.max()
wavest = wavest0.copy()
# Display wavelet
fig, axs = plt.subplots(1, 2, figsize=(20, 5))
fig.suptitle('Statistical wavelet estimate')
axs[0].plot(fwest[:nfft//2], wavest0_fft[:nfft//2], 'k')
axs[0].set_title('Frequency')
axs[1].plot(twav, wav, 'r', label='Wav true')
axs[1].plot(twav, wavest, 'k', label='Wav statistical')
axs[1].legend()
axs[1].set_title('Time');
In [66]:
fd_plus = fd_plus_sparse.copy()
d = d_sparse.copy()
In [67]:
sigmas = np.zeros(nouter)
taus = np.zeros(nouter)
scs = np.zeros((nouter, 2))
d = d.ravel()
f1_alt_inv = np.zeros(2*nt2 * nr)
tau = 0
sigma = np.inf
iouter = 0
while iouter < nouter and sigma > optimalsigma:
# 1. compute current misfit
Wavop = Convolve1D(ns*nt2, wavest, wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Wav1op = Convolve1D(ns*nt2, np.flipud(wavest), wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Mop = VStack([HStack([Wavop, -1*Wop*Rop]),
HStack([-1*Wop*Rollop*R1op, Wav1op])])*BlockDiag([Wop, Wop])
dw = np.zeros(2*nt2 * nr) if iouter == 0 else Mop * f1_alt_inv.ravel()
res = d - dw
sigma = np.linalg.norm(res.ravel())
sigmas[iouter] = sigma
# 2. update tau
lamda = - np.linalg.norm(Mop.H * res, ord=np.inf) / sigma
tau -= (sigma - optimalsigma) / lamda
taus[iouter] = tau
# 3. invert for focusing functions
f1_alt_inv = SPGL1(Mop, d, tau=tau,
x0=f1_alt_inv.ravel(),
iter_lim=50 if iouter == nouter-1 else 50,
verbosity=1)[0]
f1_alt_inv = f1_alt_inv.reshape(2*nt2, nr)
dest = Mop * f1_alt_inv.ravel()
res = d - dest
res = res.reshape(2*nt2, nr)
Mf = np.vstack((Mop * np.pad(f1_alt_inv[:nt2].ravel(), (0, nt2*nr),
mode='constant'),
Mop * np.pad(f1_alt_inv[nt2:].ravel(), (nt2*nr, 0),
mode='constant'))).T
# 4. scale focusing functions to account for amplitude underestimation
#sc = solve(Mf.T @ Mf, Mf.T @ d)
sc = np.linalg.lstsq(Mf, d)[0]
scs[iouter] = sc
print(sc)
f1_alt_inv_sc = f1_alt_inv.ravel() * np.hstack((sc[0]*np.ones(nt2*nr), sc[1]*np.ones(nt2*nr)))
f1_alt_inv_sc = f1_alt_inv_sc.reshape(2*nt2, nr)
dest_sc = Mop * f1_alt_inv_sc.ravel()
res_sc = d - dest_sc
res_sc = res_sc.reshape(2*nt2, nr)
print('res l2', np.linalg.norm(res.ravel()))
print('res_sc l2', np.linalg.norm(res_sc.ravel()))
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15, 9))
axs[0].imshow(d.reshape(2*nt2, nr), cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('data')
axs[0].axis('tight')
axs[1].imshow(res.reshape(2*nt2, nr), cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title('res')
axs[1].axis('tight')
axs[2].imshow(res_sc.reshape(2*nt2, nr), cmap='gray',
vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title('res scaled')
axs[2].axis('tight')
f1_alt_inv_tot = f1_alt_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_alt_inv_minus, f1_alt_inv_plus = f1_alt_inv_tot[:nt2].T, f1_alt_inv_tot[nt2:].T
f1_alt_inv_sc_tot = f1_alt_inv_sc + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_alt_inv_sc_minus, f1_alt_inv_sc_plus = f1_alt_inv_sc_tot[:nt2].T, f1_alt_inv_sc_tot[nt2:].T
if quickplot:
f1_alt_inv_tot_plot = np.concatenate((f1_alt_inv_tot[:(2*nt-1)], f1_alt_inv_tot[(2*nt-1):]), axis=1).T
fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 7))
ax.imshow(f1_alt_inv_tot_plot.T, cmap='gray',
vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
ax.set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
ax.axis('tight')
ax.set_ylim(1, -1);
# 5. invert for wavelet
W1op = Diagonal(w.flatten())
F1plus = VStack([Convolve1D(nt2, f1_alt_inv_sc[nt2:, ir], offset=nt-2) for ir in range(nr)])
F1minus = VStack([Convolve1D(nt2, f1_alt_inv_sc[:nt2, ir], offset=nt-1) for ir in range(nr)])
S = Flip(nt2)
Fw = VStack([W1op*F1minus, W1op*F1plus*S])
d1 = Wop*Rop*f1_alt_inv_sc_plus.T.flatten()
d2 = Wop*R1op*f1_alt_inv_sc_minus.T.flatten()
dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))
wavest = np.pad(wavest, (nt-ntwav, nt-ntwav), mode='constant')
wavest = lsqr(Fw, dw.flatten(), x0=wavest, iter_lim=20, show=False)[0]
resw = dw.flatten() - Fw @ wavest
resw = resw.reshape(2*ns, nt2)
wavest = wavest[nt-ntwav:nt+ntwav-1]
if quickplot:
plt.figure()
plt.plot(twav, wav, 'k', lw=2, label='True')
plt.plot(twav, wavest, 'r', label='Est')
plt.plot(twav, wavest0, '--b', label='Est it=0')
plt.legend()
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(5, 7))
axs[0].imshow(dw.T, cmap='gray',
vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$d_{w}$'), axs[0].set_xlabel(r'$x_R$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(resw.T, cmap='gray',
vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$r_{w}$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
iouter += 1
In [68]:
g_alt_inv = Gop*f1_alt_inv_tot.flatten()
g_alt_inv = g_alt_inv.reshape(2*nt2, ns)
f1_sparse_inv_minus, f1_sparse_inv_plus = f1_sparse_inv_tot[:nt2].T, f1_sparse_inv_tot[nt2:].T
g_alt_inv_minus, g_alt_inv_plus = -g_alt_inv[:nt2].T, np.fliplr(g_alt_inv[nt2:].T)
g_alt_inv_tot = g_alt_inv_minus + g_alt_inv_plus
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_alt_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_alt_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_alt_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_alt_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_alt_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip,
extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_alt_inv_tot[nr//2, nt-1:]/g_alt_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
In [69]:
fig, axs = plt.subplots(1, 3, figsize=(15, 4))
axs[0].plot(taus, 'k')
axs[0].set_title('tau')
axs[1].plot(scs[:, 0], 'k', label='1st term')
axs[1].plot(scs[:, 1], '--k', label='2nd term')
axs[1].set_title('sc')
axs[2].plot(sigmas, 'k')
axs[2].axhline(optimalsigma, color='k', lw=5)
axs[2].set_title('sigma');
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
ax.plot(taus, sigmas, 'k')
ax.title('tau vs sigma')
iouter=0 tau=0
d = d.ravel() f1_alt_inv = np.zeros(2nt2 nr)
1. compute current misfit¶
Wavop = Convolve1D(nsnt2, wavest, wav_c, dims=(nt2, ns), dir=0, dtype='float64') Wav1op = Convolve1D(nsnt2, np.flipud(wavest), wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Mop = VStack([HStack([Wavop, -1WopRop]), HStack([-1WopRollopR1op, Wav1op])])BlockDiag([Wop, Wop])
dw = np.zeros(2nt2 nr) if iouter == 0 else Mop * f1_alt_inv.ravel() res = d - dw sigma = np.linalg.norm(res.ravel())
2. update tau¶
lamda = - np.linalg.norm(Mop.H * res, ord=np.inf) / sigma tau -= (sigma - optimalsigma) / lamda
3. invert for focusing functions¶
f1_alt_inv = SPGL1(Mop, d, tau=tau, iter_lim=10, verbosity=1)[0] # x0=f1_alt_inv.ravel(), f1_alt_inv = f1_alt_inv.reshape(2nt2, nr) dest = Mop f1_alt_inv.ravel()
Mf = np.vstack((Mop np.pad(f1_alt_inv[:nt2].ravel(), (0, nt2nr), mode='constant'), Mop np.pad(f1_alt_inv[nt2:].ravel(), (nt2nr, 0), mode='constant'))).T
4. scale focusing functions to account for amplitude underestimation¶
sc = solve(Mf.T @ Mf, Mf.T @ d)¶
sc = np.linalg.lstsq(Mf, d)[0] print(sc) f1_alt_inv = f1_alt_inv.ravel() np.hstack((sc[0]np.ones(nt2nr), sc[1]np.ones(nt2nr))) f1_alt_inv = f1_alt_inv.reshape(2nt2, nr)
dest_sc = Mop * f1_alt_inv.ravel()
ffig, axs = plt.subplots(1, 3, sharey=True, figsize=(15, 9)) axs[0].imshow(d.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[0].axis('tight') axs[1].imshow(dest.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[1].axis('tight') axs[2].imshow(dest_sc.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[2].axis('tight')
f1_alt_inv_tot = f1_alt_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus)) f1_alt_inv_minus, f1_alt_inv_plus = f1_alt_inv_tot[:nt2].T, f1_alt_inv_tot[nt2:].T
f1_alt_inv_tot_plot = np.concatenate((f1_alt_inv_tot[:(2nt-1)], f1_alt_inv_tot[(2nt-1):]), axis=1).T fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 7)) ax.imshow(f1_alt_inv_tot_plot.T, cmap='gray', vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1])) ax.settitle(r'$f{inv}$'), axs[0].set_xlabel(r'$x_R$') ax.axis('tight') ax.set_ylim(1, -1);
5. invert for wavelet¶
W1op = Diagonal(w.flatten()) F1plus = VStack([Convolve1D(nt2, f1_alt_inv[nt2:, ir], offset=nt-2) for ir in range(nr)]) F1minus = VStack([Convolve1D(nt2, f1_alt_inv[:nt2, ir], offset=nt-1) for ir in range(nr)]) S = Flip(nt2) Fw = VStack([W1opF1minus, W1opF1plus*S])
d1 = WopRopf1_alt_inv_plus.T.flatten() d2 = WopR1opf1_alt_inv_minus.T.flatten() dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))
wavest = np.pad(wavest, (nt-ntwav, nt-ntwav), mode='constant') wavest = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0] # x0=wavest wavest = wavest[nt-ntwav:nt+ntwav-1]
if quickplot: plt.figure() plt.plot(twav, wav, 'k', lw=2) plt.plot(twav, wavest, '--r');
6. remove scaling from focusing functions¶
f1_alt_inv = f1_alt_inv.ravel() np.hstack((sc[0]np.ones(nt2nr), sc[1]np.ones(nt2nr))) f1_alt_inv = f1_alt_inv.reshape(2nt2, nr)
In [ ]: